home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / tmtdemo_r < prev    next >
Encoding:
Text File  |  1995-01-17  |  4.5 KB  |  179 lines
  1. //
  2. // TESTDEMO       Demonstration of Test Matrix Toolbox.
  3. //                N. J. Higham.
  4.  
  5. // This file is a translation of tmtdemo.m from version 2.0 of
  6. // "The Test Matrix Toolbox for Matlab", described in Numerical
  7. // Analysis Report No. 237, December 1993, by N. J. Higham.
  8.  
  9. "The Test Matrix Toolbox contains test matrices, visualization routines,"
  10. "and other miscellaneous routines.\n"
  11.  
  12. "The version of the toolbox is RLaB Port of v2.0\n"
  13.  
  14. "For this demonstration you will need to view both the command window"
  15. "and one figure window."
  16. "This demonstration emphasises graphics and shows only"
  17. "some of the features of the toolbox.\n"
  18.  
  19. "First, lets load all of the Rfiles in the testmatrix directory"
  20. "Note: you must be executing this file from the testmatrix directory.\n"
  21.  
  22. rfile loadtm
  23.  
  24. "The FV command plots the boundary of the field of values of a matrix"
  25. "(the set of all Rayleigh quotients) and plots the eigenvalues as"
  26. "crosses (`x').  Here are some examples:\n"
  27.  
  28. "The Grcar matrix is a Toeplitz matrix of the following form:"
  29.  
  30. grcar(5)
  31.  
  32. "Here is the field of values of the 10-by-10 Grcar matrix:"
  33.  
  34. ptitle ("fv (grcar (10))");
  35. fv (grcar (10));
  36. pause ("");
  37.  
  38. "Next, we form a random orthogonal matrix and look at its field of values."
  39. "The boundary is the convex hull of the eigenvalues since A is normal.\n"
  40.  
  41. ptitle("randsvd (10, 1)");
  42. A = randsvd (10, 1);
  43. fv (A);
  44. pause ();
  45.  
  46. "The RANDSVD commands generates random matrices with pre-assigned"
  47. "condition number, and various singular value distributions:\n"
  48.  
  49. A = randsvd(6, 1e6, 3);  // Exponential distribution.
  50.  
  51. svd(A).sigma.'
  52. 1/rcond (A)
  53. pause ();
  54.  
  55. "The PS command plots an approximation to a pseudospectrum of A,"
  56. "which is the set of complex numbers that are eigenvalues of some"
  57. "perturbed matrix A + E, with the norm of E at most epsilon"
  58. "(default: epsilon = 1E-3)."
  59. "The eigenvalues of A are plotted as crosses (`x')."
  60. "Here are some interesting PS plots.\n"
  61.  
  62. "First, we use the KAHAN matrix, a triangular matrix made up of sines and"
  63. "cosines.  Here is an approximate pseudospectrum of the 10-by-10 matrix:\n"
  64.  
  65. ptitle ("ps(kahan(10))");
  66. ps (kahan (10), 25);
  67. pause ();
  68.  
  69. "The TRIW matrix is upper triangular, made up of 1s and -1s:\n"
  70.  
  71. triw (4)
  72.  
  73. "The following matrix has a pseudospectrum in the form of a limacon.\n"
  74.  
  75. ptitle ();
  76. n = 25; 
  77. A = triw(n,1,2) - eye(n,n);
  78. sub (A, 6)               // Leading principal 6-by-6 submatrix of A.
  79. ps (A);
  80. pause ();
  81.  
  82. "Here is the 8-by-8 Frank matrix.\n"
  83. A = frank(8)
  84.  
  85. "We can get a visual representation of the matrix using the SEE"
  86. "command, which produces subplots with the following layout:"
  87. "    |---------------------------------|"
  88. "    | MESH(A)            MESH(INV(A)) |"
  89. "    | SEMILOGY(SVD(A))   FV(A)        |"
  90. "    |---------------------------------|"
  91. "where FV is the field of values.\n"
  92.  
  93. pclose ();
  94. pstart (2,2);
  95. see(A);
  96. pause ();
  97.  
  98. "The Frank matrix is well-known for having ill-conditioned eigenvalues."
  99. "Here are the eigenvalues (in column 1) together with the corresponding"
  100. "eigenvalue condition numbers (in column 2):\n"
  101.  
  102. es = eigsens (A);
  103. [es.val.', es.sens, 1./es.val.']
  104.  
  105. "In the last column are shown the reciprocals of the eigenvalues."
  106. "Notice that if LAMBDA is an eigenvalue, so is 1/LAMBDA!\n"
  107.  
  108. pause ();
  109.  
  110. "MAGIC function produces magic squares:\n"
  111.  
  112. A = magic(5)
  113.  
  114. "Using the toolbox routine PNORM we can estimate the matrix p-norm"
  115. "for any value of p.\n"
  116.  
  117. [pnorm(A,1), pnorm(A,1.5), pnorm(A,2), pnorm(A,pi), pnorm(A,inf())]
  118.  
  119. "As this example suggests, the p-norm of a magic square is"
  120. "constant for all p!\n"
  121.  
  122. pause ();
  123.  
  124. "GERSH plots Gershgorin disks.  Here are some interesting examples.\n"
  125.  
  126. ptitle("gersh(lesp(12))");
  127. gersh(lesp(12));
  128. pause ();
  129.  
  130.  
  131. ptitle("gersh(hanowa(10))");
  132. gersh(hanowa(10));
  133. pause ();
  134.  
  135. ptitle("gersh(ipjfact(6,1))");
  136. gersh(ipjfact(6,1).A);
  137. pause ();
  138.  
  139. ptitle("gersh(smoke(16,1))");
  140. gersh(smoke(16,1));
  141. pause ();
  142. pclose();
  143.  
  144.  
  145. "A Hadamard matrix has elements 1 or -1 and mutually orthogonal rows:"
  146.  
  147. hadamard(4)
  148.  
  149. "A CONTOUR plot of this matrix is interesting:"
  150.  
  151. pstart (1,1);
  152. plcont(<< x=1:16; y=1:16; z=hadamard(16) >>);
  153. pause ();
  154.  
  155. "GFPP generates matrices for which Gaussian elimination with partial"
  156. "pivoting produces a large growth factor."
  157.  
  158. gfpp(6)
  159. pause ();
  160.  
  161. "Let's find the growth factor for partial pivoting and complete pivoting"
  162. "for a bigger matrix:"
  163.  
  164. A = gfpp(20);
  165.  
  166. "Partial pivoting."
  167. LU = lu(A);            
  168. max(max(abs(LU.u))) / max(max(abs(A)))
  169.  
  170. "Complete pivoting using toolbox routine GECP."
  171. LUPQrho = gecp(A);     
  172. LUPQrho.rho
  173.  
  174. "As expected, complete pivoting does not produce large growth here."
  175.  
  176. pause ();
  177. clearall ();
  178. format ();
  179.